On the gluon spectrum in the glasma 
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I. INTRODUCTION 

In the Color Glass Condensate (CGC) framework (for 
reviews see [ll, 0), particle multiplicities are dominated 
by classical gluon dynamics. The spectra of gluons pro- 
duced in the collision of two heavy nuclei at high energy 
is determined [J-IH completely by the classical glasma [q] 
gauge field A^^ that is the solution of the classical Yang- 
Mills (CYM) equations of motion: 



[D^,Fn = r, 



(1) 



d, 



fi — igAfj, is the covariant derivative and 



where D^ 

Fij,„ = dfj,A^ — duAfi — ig[A^, A^,]. The current J'', which 
must be covariantly conserved ([D^, J''] =0), describes 
the colliding nuclei. At lowest order in the color densities 
of the nuclei it reads 



J^ 



p^(a;+,XT), J- = P^(a;-,XT), J' = 0, (2) 



with x^ — {t ± z)/^/2. This current reflects the kine- 
matics of two nuclei A and B moving in the —z and 
+z direction, respectively, at almost the speed of light. 
With a specific gauge choice, for instance with the light- 
cone gauge or the Fock-Schwinger gauge, the conserved 
current can be fully determined, allowing us to solve 
the CYM equations on the light-cone. The problem of 
gluon production reduces then to an initial value prob- 
lem of solving the CYM equations [D^j„F^''] = in the 
forward light-cone, x > 0. The initial condition for 
the gauge field can be determined explicitly by solv- 
ing the CYM near the light-cone, i.e., at proper time 



v^ 



= 0+ 



An often used approximation for computing gluon pro- 
duction in heavy ion collisions is to express the spectrum 
as a convolution of the unintegrated gluon distribution 
functions, leading to the so-called /cr-factorization, which 
holds exactly in the dilute "pp" (proton-proton) and the 
dilute-dense "pA" (proton-nucleus) limits, that is, when 
one or both of the color sources are weak |7Hl0l| . In these 
cases the CYM equations can be linearized and solved 



explicitly^ . 

In a recent work [12| the spectrum of gluons has been 
obtained by solving the Yang-Mills equations near the 
light-cone and treating the evolution in time linearly. It 
has been shown that this procedure yields the correct 
"pA" spectrum after expanding the solution at first order 
in the weak proton sources, provided one chooses a gauge 
that removes the transverse pure gauge component of 
the fields of the nuclei when they are not interacting. 
In the present paper, we investigate the accuracy of this 
approach in the nucleus-nucleus case. More precisely, the 
following questions will be addressed: 

• What is the role of the nonlinearities in the initial 
condition of the CYM calculations and in the time 
evolution for r > 0? How well can one approximate 
the full solution of the nonlinear equations by one 
where the initial condition is treated exactly but 
the time evolution is assumed linear? 

• How good an approximation is a /cr-factorized ex- 
pression for the gluon spectrum in the glasma? 

We will argue in the following that a specific gauge, 
namely the Fock-Schwinger gauge A^- — combined with 
a two dimensional transverse Coulomb gauge condition 
diAi = at T = 0"*", allows us to minimize the magni- 
tude of the initial gauge fields and therefore to reduce the 
"final state interactions," i.e., the non- linear dynamics 
inside the forward-light-cone. We shall then, in Sec. |Vl 
discuss the range of validity of fcr-factorization. 

In order to answer the questions above, we shall com- 
pare approximate solutions with the full numerical solu- 
tions of the CYM equations, obtained using the method 
developed in Refs. |13l - [l7j . As we shall discuss in Sec.lTll 
the initial conditions for the solution are given in terms 



mject 



For a conjectured analytical solution in the fully nonlinear case 
see Ref. 



of Wilson lines that, in the CGC, are random SU(3) ma- 
trices, correlated in the transverse plane over distances 
of order l/Qs- For our numerical studies we shall gen- 
erate these Wilson lines from an implementation of the 
McLerran-Ve nugo palan (MV) model [j-Q (for more de- 
tail see Refs. [la, Ss U^), where the color densities are 
treated as random variables with a Gaussian distribution 
of variance 

g^t,\x+)5^'S{x+~y+)5i^T-yT), (3) 

for nucleus A and similarly for nucleus B. Here, g^jjp'^x^) 
stands for the three-dimensional (squared) color charge 
density. After the averaging over the color densities, the 
observable will be only function of the two dimensional 
density of (squared) color charge 



2 2 

^ A*. = 



+ 00 



dx+gV^(-^+) 



(4) 



We denote the value integrated over the longitudinal co- 
ordinate X by /i^ without an explicit dependence on 
X . The saturation scale Qs is is proportional to (and 
numerically of the same order of magnitude as) g^lJ-^- 
The relevant values are of order g^/i^ — 1 — 2 GcV for 
RHIC and 2 - 3 GeV for the LHC. Note that in the 
convention used here ii\ is a number density of charge 
carriers and ff^/x^ the squared charge density. In the sat- 
uration regime ff^M^ ~ l/o^s so that the saturation scale 
Qs ~ ff^/^A is parametrically independent of the coupling. 
For a discussion of the relation between the MV model 
parameter g^fjL and the saturation scale Qs as measured in 
DIS experiments we refer the reader to Refs. [l^,[l^. We 
emphasize that the main purpose of this paper is not to 
compare directly to experimental data, but to compare 
different approximations for the gluon spectrum in the 
glasma corresponding to the same distribution of color 
charges (or more properly Wilson lines, as we shall see 
shortly). Our main results would therefore be indepen- 
dent of the details of the color charge density distribution, 
and equally valid for one that would match more closely, 
e.g., the solution of the JIMWLK/BK evolution equa- 
tions. Note however, that although we do not attempt 
to choose values for the parameters of our calculation to 
give an accurate best estimate for RHIC or the LHC, we 
shall stay in the phenomenologically relevant range. 



II. LINEARIZED APPROXIMATION OF 
YANG-MILLS EQUATIONS 

In this section we describe our linearized approxima- 
tion for solving the GYM equations. 

We work in the the Fock-Schwinger (FS) gauge tAt = 



-A- 



-A+ =0. InRef. 



the discussion was formu- 



lated in terms of the light cone (LG) gauge. Appendix G 
of Ref. [l2] shows how both gauge conditions lead to the 



same functional form of the spectrum in the linearized 
approximation, Eq. (P5)) of the present paper. The re- 
maining longitudinal component of the gauge field is most 
naturally parametrized by the component A^ orthogonal 
to At, namely as A^ — ix^A^, where "H — \ \n{x'^ /x~) 
is the spacetime rapidity. 

The initial condition for the field on the light-cone is 
given by i2a^l] 



4M — 4' 



■Al 



A" 



«5r 



T=o+ - 2 ^ A' £ 



(5) 



where the LC gauge fields for the individual nuclei are 
given by 



A\ 



—U'^d'U, A' 
ig 



-V^d'V. 



W 



(6) 



The fundamental degrees of freedom characterizing the 
GGC wavefunctions of the individual nuclei are the Wil- 
son lines C/(xt) and ^(xt). 



C/(xt) = P+ exp 
and 

V{xt) = V- exp 



*5 



/ dz+— 2-p^(z+,xt) 



ig / dz ■=2Pb{z ,xt) 



(7) 



(8) 



where p^(x+,xt) and p^(x~,xt) are the color charge 
densities of the nuclei. In the numerical implementation 
of the MV model these Wilson lines are constructed as 



t/(xT) 



Ny 

J|exp 

fc=i 



1 



-«ff: 



■pn^T] 



(9) 



where the color charges are Gaussian variables with the 
variance 

{pU^T)pHyT)) = s^'6'^'s'iKT - yr)^. (lo) 

The indices fc, / = 1 . . . Ny represent a discretization of 
the longitudinal direction into Ny small steps; the contin- 
uum limit corresponding to Eqs. (O and (|S]) is achieved 
for Ny — >■ oo at constant g^ [i. Some kind of infrared 
regulator is needed in order to invert the Laplacian oper- 
ator V^. One possibility is to replace V^ in Eq. © by 
Vy -f m^ , with a parameter m chosen so that m <^ Qs- 
Another possibility is to set the ky = Or-mode of p to 
zero, i.e. to impose total color neutrality on the source. 
This corresponds to an infrared cutoff given by the size 
of the system; this is the procedure we use in what fol- 
lows in the cases where we take m — 0. As noted before, 
our purpose is to compare different methods to obtain 
the gluon spectrum from the same distribution of Wilson 
lines. Therefore the important comparison in this context 
is between different approximations for the same values 



of Ny, m and g^/i. We shall not carry out a a systematic 
study of the dependence on these parameters separately 
(see Ref. [l^ for a more detailed analysis), except for the 
specific case of the dependence on Ny in Sec. IIVI 

Note that when the fields are boost invariant, the 
gauge condition At- ~ Q leads to d+A^ + 9_y4^ = 0. 
This gauge, as well as the LC-gauge, does not fix com- 
pletely the gauge field, and we may exploit the remain- 
ing gauge freedom in order to simplify the calculation. 
Imposing the restriction that we want to stay within the 
Fock-Schwinger gauge A^. —Q forbids r-dependent gauge 
transformations. We also want to preserve the explicit 
boost invariance of the field configurations and therefore 
we do not want to perform ?7-dependent transformations. 
This leaves us the freedom of gauge changes that depend 
(in the region r > 0) only on the transverse coordinates. 
We denote the transformed field by A: 



A^ = nA^" n'^ - —^d^'nK 
19 

The initial conditions in the new gauge are 



(11) 



A\ 



-0- = |^[Al,A;]f7t 



=0+ 



= n(A\ + Ai)n^ - —nd'n\ (12) 

ig 



where A\ and A^ are given by Eq. 



As we shall see. 



we can choose il so as to reduce final state interactions 
and treat the evolution of the fields after the collision in 
a linear approximation. As recalled earlier, a motivation 
for this strategy is the fact that, in a suitable gauge, 
treating the final state dynamics to lowest order in the 
field gives the exact solution in the case where one of the 
projectile is dilute and the other one dense (the "pA" 
case). 

We shall return to the determination of il later, and 
first review the linearized solution, following Ref. [12l |. 
The linearized equations of motion, for r > 0, are 



DA' 
DA^ 



-d'd'A, 
-d^d^A. 



(13) 
(14) 



These are to be solved with the initial conditions given 
in Eq. p^ . The second equation (ITil) . combined with 
the gauge condition and boost invariance, leads to 

D [d+A- + d-A+] = 2d+d-d^A = 0, (15) 

which states that the divergence of the transverse field 
is conserved. This allows us to rewrite the Yang-Mills 
equations in a form that makes the boundary conditions 
explicit. Following the same steps as in Ref. [l2|, we 
obtain 



DA' 

DA+ 
DA- 



25{x+)Six-) Al^^^ 

^e{x+)eix-)d'd^ A\^ 
s{x-)e{x+) A'^U,^ , 

~Six+)9{x~) A\^o+ 



=0+ ' 



(16) 

(17) 
(18) 



Fourier transforming these equations we get 



k'Aik) = 2[6'^ 



2k+k' 



- -A'O^t) 



T=0 + 



-fcM+(fc) = -— ^"(MU^, 



^k^A-{k) 



A^k^ 



It=0+ 



,(19) 
(20) 
(21) 



The gluon spectrum is then computed with the help of 
the reduction formula: 



An'^E 



dN 

d3k 



Y^\M> 



(22) 



with the production amplitude Ai\ — lim k^e'^A^{k) 
for a gluon of polarization A. By using the completeness 

E^( 



relation J2^n{f'u)* = ^.9^"^' ^^^ finally gets [12J : 



dN 



1 



1 



dyd2kT (27r)2 7rk2, 



|kT X A{kT)r 

+ [^"(kr) 



T = 0+ 



(23) 



where the cross product stands for k^ x ^ = e'^k'^A^ . 

Let us now return to the choice of the gauge transfor- 
mation n in Eq. (fTTT) . The choice introduced in Ref. pj[ 
(see Appendix C of [1^ for the explicit derivation) was to 
take either O = UV or Q = VU. It was shown in Ref. [IJl 
that this choice reproduces correctly both the "pA" and 
"Ap" limits (because the gauge is not symmetric in the 
two nuclei one has to study the two cases separately). 
In general, any choice that would gauge away the pure 
transverse fields of the nuclei, A^ and A^ , before the 
collision, leads to the correct answer for "pA". Thus we 
look for fl that verify the following boundary conditions: 
when x^ < 0, i.e., in the absence of nucleus B, fi = [/ and 
when a;+ < 0, i.e., in the absence of nucleus A, J7 = V. 

With the explicit expression Jl = VU one can directly 
perform the gauge rotations of the initial conditions in 
Eq. (dH) to get^ 

■^1=0+ = V{UA^U^^A^)V^ (24) 

A^l^^Q., = V{d'UA^U^ + UA^d'U'^)V'<. (25) 

The expression for the spectrum that results from in- 
serting Eqs. dM]) and ^5^ in Eq. ^^ was obtained in 
Ref. [ijj. 

We now look for a gauge where the linearized evolution 
for T > works best, that is, we look for a gauge that 
minimizes the value of the gauge potential A. With our 
restriction of boost invariance and the Fock-Schwinger 



2 These are Eqs. (2.73) and (2.74) of Ref. [TJ] written in the 
fundamental representation. 



condition Ar = a natural way to do this is to minimize 
the transverse components of A. They are the ones in 
which the large unphysical pure gauge contributions show 
up. as can be explicitly seen in the computation of Refs. 
[2Ct |22| . A convenient prescription for this is to minimize 
the functional / d'^xrl At {xt,t — 0^)^; this is achieved 
by the 2-dimensional Coulomb gauge d^A^ = 0. The 
Coulomb gauge condition is conserved by the linearized 
equations of motion (see Eq. p5|) ). so it is equivalent to 
impose it at r = 0+ or at a larger r. The full equations 
of motion, on the other hand, do not conserve the gauge 
condition. Thus in the numerical computation the gauge 
condition is imposed at the end of the evolution, when 
the gluon spectrum is calculated. 

In order to impose the condition d^A^ — wc must 
then find the gauge transformation il as a solution of the 
equation 



d^ n{A^^+Ai) 



Ai)n^ ~—nd'n^ 



w 



(26) 



When x^ < Q, VL = U and when x+ < Q, Vl = V as 
required to reproduce the "pA" spectrum. In this case 
the pure transverse initial fields get rotated leading to 
pure longitudinal fields, i.e.. 



A-" 



1 



-p^ {x , xt), a = A' = 0, for x+ < 0, 



A^ = ^P^{x+,x.t), A'^ =A' =0, ioTx < 0. 

(27) 

In the strong field case we are not able to solve Eq. ([25)1 
for fl analytically (or even to show formally that a unique 
solution exists). Finding the required gauge transforma- 
tion numerically is, however, not excessively difficult and 
has been a standard part of the numerical CYM compu- 
tations of the gluon multiplicity T4l-4i3!|. Indeed evalu- 
ating the approximate formula Eq. (j23p in the Coulomb 
gauge defined by Eq. ([^5]) is significantly less demanding 
than solving the full time-dependence of the Yang-Mills 
equations. 

In the following sections we shall compare the approx- 
imate analytic solution Eq. (^5)) in the two gauges dis- 
cussed above, ^ — VU, and the Coulomb gauge choice, 
to the full numerical solution of the Yang-Mills equations. 
In addition to giving physical insight into the (admittedly 
gauge dependent) question of the importance of initial 
(meaning r = 0) and final (r '^ 1/Qs in this context) 
state interactions in the glasma this may also constitute 
a starting point for further analytical studies. 



III. LINEAR AND NONLINEAR FINAL STATE 
DYNAMICS 

We start by comparing the results of the linearized evo- 
lution in the il = VU case and the Coulomb case in the 



— g |ix = 12 

- - X = 0, Coul. 
x = 0,VU 




FIG. I: il = VU gauge gauge and Coulomb gauge results 
compared to the full CYM result in the weak field regime 
with gV = 0.2 GeV, m = 0.1 GeV, Ny = 20. The legend 
"r = 0" refers to the spectrum being evaluated using Eq. pS)) 
which expresses the multiplicity in terms of the fields at r = 0. 
The full CYM result is evaluated at r = 12/gV- 




FIG. 2: Q, = VU gauge and Coulomb gauge results compared 
to the full CYM result in the "pA" case g'^fi^ = 1.25 GeV, 
g fj,p — 0.008 GeV. There are two n = VU gauge curves, 
corresponding to the cases where either the source A or B 
is taken to be weak, i.e. V^Up: source A is the proton, and 
VpU^: source B is the proton. It can be seen that in this 
limit of small g^ ^ip they are equivalent, but the Vg Up approx- 
imation approaches the dilute limit more slowly. The CYM 
result is evaluated at r = 60/(?^/i^. Ny — 20 and m, = 0. 



dilute linfit for both sources. Figure [T] shows the gluon 
spectrum for .g^/i = 0.2 GeV, Ny = 20 and m = 0.1 GeV. 
This situation is dilute enough to find a good agreement; 
the result mainly serves as a check of the normalization 
in our numerical computation. The uneven structure at 
small kr in the full CYM calculation is an oscillation 
caused by the fact that the evolution is stopped at a finite 
time r; see Appendix ID] for a more detailed discussion. 



g |ix = 12 
X = 0, Coul. 
x = 0,VU 




FIG. 3: Light-cone gauge and Coulomb gauge results com- 
pared to the full GYM result in the saturated strong field 
regime with g'^/i = 2 GeV, m = 0.1 GeV A''^ = 20. 



A slightly less trivial check is provided by the "pA" 
case. We write this in quotation marks because what 
we are really computing is the dilute-dense limit where 
the smaller one of the two saturation scales is taken to 
be much smaller than the other, but also the weaker 
source is allowed to fill the whole transverse plane (of 
area S± — TrR\) so that we need not worry about ef- 
fects of the edge of the nucleus. Here also several an- 
alytic calculations in the Coulomb [22], covariant 0,0] 
and LC [9, 12] gauges have shown that the result can 
be written in a fcy-factorized form, in which one needs 
to account for nonlinear interactions only in the initial 
condition. A practical question is then of course how 
weak the "proton" source has to be to be in this limit. 
Our numerical comparison is shown in Fig. [51 where the 
Coulomb and fl — VU gauge approximations are seen 
to agree with the full numerical computation. Note that 
because the fi = VU gauge prescription, Eqs. ([24]) and 
(1551) . is asymmetric in the nuclei A and B, one has to 
look at the "pA" and "Ap" limits of the approximation 
separately. In Fig. [2] these are denoted by Vg Up (source 
A is the proton) and VpUJ^ (source B is the proton). In- 
specting the figure closely one can see that the VgUp 
curve is slightly further away from the full CYM result. 
This signals the fact that the limit of nucleus B being di- 
lute converges towards the linear regime faster than the 
limit of A being dilute, which we have verified for less 
asymmetric values of g^fi^ and g^^ip- In other words, 
for the "marginal pA" case in which the proton source 
is weak but not infinitesimal, the gauge transformation 
Vt = VpUj^ gives a slightly better approximation to the 
full result than f2 — V^Up, where the Wilson line of the 
proton is denoted with subscript p ■^. 



Let us then turn to the same comparison for the 
nucleus- nucleus case when both sources are strong. The 
comparison between the linearized prescriptions and the 
full result is shown in Fig. [31 We can see that the lin- 
earized Coulomb gauge result is still relatively close to the 
full numerical solution, but the one in VU-gaiige starts 
to deviate from it significantly. After our discussion in 
Sec. [n] this behavior is easy to interpret. The Coulomb 
gauge condition removes, practically by construction, the 
unphysical pure gauge component from the field. A lin- 
earized approximation works best in a gauge that min- 
imizes the value of the gauge potential. The T^[/-gauge 
condition removes the pure gauge part of the field in the 
"pA" case, but not in the "AA" one, because the gauge 
transformation required to do this is more complicated. 
Let us comment in more details the comparison between 
the linearized Coulomb gauge result with the full CYM 
result. We distinguish typically three regions in the spec- 
trum: 

• The high momentum range, kx/g^fJ- > 3. The two 
curves merge together, this is expected since both 
tend to the same limit ^ l/fc|, at high momentum. 

• The intermediate momentum region, 0.2 < 
kx/g^fj. < 3 , where the two curves follow the same 
trend and are peaked at kx/g^^J- — 1.2. In this re- 
gion, our analytical formula overestimates the full 
CYM result by less than 15 %. 

• The low momentum region, kr/g^J^ < 0.2. Here, 
we see that the two curves diverge. The full result 
tends to whereas our approximate formula tends 
to a constant. 

Note that the spectrum is multiplied by k^ in our figures, 
which focuses attention on the larger fcr-part and min- 
imizes the difference between the exact calculation and 
the approximate ones. But the difference between the 
Coulomb gauge result and the full CYM result for the 
integrated multiplicity remains large although the spec- 
tra are close in shape. The difference in the integrated 
multiplicity comes mostly from the infrared part of the 
spectrum. Indeed since the Coulomb-gauge spectrum is 
logarithmically IR divergent it does not give a reliable 
estimate for the integrated multiplicity. One may ar- 
gue that, due to the uncertainty in the CYM calculation 
for kx ^ 1/t the very small fc^-part of the spectrum 
in Fig. [31 is is not very reliable either. Earlier numeri- 
cal studies jl4)-tlq have shown, however, that the total 



To trace the origin of this small effect, look at the analytical ex- 



pressions for the "pA" and "Ap" limits, Eqs. (3.91) and (3.96), 
in Ref. [1^ . In Eq. (3.91) px is parametrically large, ^^ QsA 
and ky small ~ Qsp- The longitudinal component of the gauge 
field a' is proportional to the large momentum p'. In Eq. (3.96) 
ky is large and p^ small, and the longitudinal component of 5* 
is proportional to the small momentum p*. Since the "unphys- 
ical" longitudinal polarization is smaller in (3.96), it is a better 
approximation. 
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FIG. 4: Nuclear modification factor of tlie gluon spectrum 
in the saturated asymmetric "pA" case ^^ — 2 GeV, g^ ^J.p ~ 



0.32 GeV, m = 0.1 GeV iV^ = 500. 



FIG. 5: Nuclear modification factor Raa, normalized with 
(g^fj.)'^ of the gluon spectrum for g^/i^ — 2 GeV, g^fJ-p = 
0.32 GeV and Ny = 500. 



gluon multiplicity is infrared finite. This IR-finitcncss 
of the CYM result, presumably due to screening effects 
at small transverse momenta [IJ, [2j, [2J] , cannot be cap- 
tured by a linearized treatment of the evolution for r > 0. 



IV. NUCLEAR MODIFICATION FACTOR 

We have shown in the previous section that an ini- 
tial condition including all orders in the classical field in 
the Coulomb gauge followed by a linearized solution to 
the equations of motion gives a good approximation to 
the full solution except at very small momenta. Let us 
now further illustrate this point by computing the nuclear 
modification factors of the gluon spectrum. These are de- 
fined as the ratio of the spectrum of produced gluons to 
the corresponding spectrum in the dilute "pp" case, nor- 
malized by the appropriate geometrical factor, the num- 
ber of binary collisions A^coU, to get a quantity that should 
approach 1 at large transverse momenta. In order to 
avoid complications with edge effects and Coulomb tails 
of the gauge field extending outside the nucleus |25l - [27| 
we compute the gluon spectrum in collisions of two ob- 
jects of the same size (filling the whole transverse lattice), 
but with different saturation scales representing a proton 
or a nucleus. We then compare the obtained gluon spec- 
tra per unit transverse area {S±) and normalize these 
using the expected large kx behavior of the gluon spec- 
trum as ^ {g'^^^)'^{g'^fig)^/k^ to get a quantity that 
approaches unity at large momenta. Thus for a collision 
between two generic nuclei A and B our definition of the 
nuclear modification factor is 



Rab 



Ai„ 



2 ,,2 



n^' 



dNAB 



dyd^kT d^Si 



AN; 



pp 



dyd^hr d^Si 



(28) 



The experimentally measured quantity is naturally inte- 
grated over the transverse area of the collision system. 
In the case of a symmetric collision Eq. (1^51) reduces to 
Raa and for a pA collision to RpA- In the MV model the 
physical interpretation of the saturation scale is straight- 
forward as resulting from an incoherent sum of the color 
charges of high- a; partons. This leads to the scaling of the 
color charge density as /i^ — A^^^^'i, and the normaliza- 
tion factor in Raa reduces to A^^^ which is simply the 
number of binary collisions per unit area dA^coii/d^S'j,. 
In the proton-nucleus case the ratio of the spectra in 
Eq. (j28p is again normalized by the number of collisions 
A^coU, which now scales like A^^^ . 

Our numerical results for RpA and Raa are shown 
in Figs. [3] and Fig. [S] respectively. They both exhibit 
a Cronin enhancement, peaked respectively at about 
two and three times the saturation scale g^/i^- The 
Cronin peak has been indeed observed in "pA" collisions, 
whereas in the "AA" case a large suppression of a factor 
5 has been observed at RHIC, commonly understood in 
terms of energy loss in the hot and dense medium pro- 
duced in the collision. The CYM calculation does not 
account for such a final sate effect. 

There is an additional remark we must make concern- 
ing the numerical calculation. In order for the nuclear 
modification factor to approach one for large transverse 
momenta, the unintegrated gluon distribution (or Wilson 
line correlator A/'(kT), see Eq. (|5D|) below) must approach 
the asymptotic behavior ~ (5^/i)/k|i with a constant of 
proportionality that is independent of any infrared scale 
in the problem. We have numerically found (see Ref. [I^l 
for more details) that to achieve this one must discretize 
the longitudinal coordinate in constructing the Wilson 
line on a very fine grid. This means that to approach the 
right large ky limit in nuclear modification factor one 
must use a much larger value of Ny than was needed for 
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Scaled fundamental representation Wilson line cor- 



relators {k^/{g^fiY)Af(]iT/g^fJ,) for g^^I^ = 2 GeV, ff>p = 



0.32 GeV, m = 0.1 GeV. 



Results are shown for both Ny = 20 



and Ny — 500; the latter ones are obtained from the configura- 
tions used for Figs. |4] and [5] The scaling A/'(kT) ~ {g'^fif/kT 
is achieved only for the larger value of Ny. When the hori- 
zontal axis is scaled with the position of the peak (i.e. Qs) 
instead of g^fj,, the difference between the different values of 



the determination of Qs in terms of g'^fi n§\ or for the 
single or double inclusive gluon spectra in the bulk region 
of momenta around Qs [3 ■ The slow convergence in the 
Ny — >■ CX3 limit is demonstrated in Fig. [5] for J\f in the 
fundamental representation. 

As we already saw in the previous section, in the dilute- 
dense limit the approximation of linearized final state 
evolution gives the correct result. For the nuclear modi- 
fication ratio RpA this is demonstrated in Fig. 21 where 
this ratio is plotted using both the full CYM result and 
the Coulomb gauge approximation. A slight difference 
between the full result and the Coulomb gauge approx- 
imation for Raa is seen in Fig. [SI the Coulomb gauge 
approximation leading to a slight overestimate. 



V. THE fcT-FACTORIZED APPROXIMATION 



As is discussed in detail in Ref. [l2|, there is no valid 
/cT-fa-ctorized expression for the gluon multiplicity in the 
fully nonlinear case of AA-collisions, in contrast to the 
case of the dilute "pA" limit. Different /cT-factorized 
approximations for nucleus-nucleus collisions have never- 
theless been extensively used in the literature and it is 
therefore instructive to check how good these approxima- 
tions are. 

Our starting point is the following /cr-factorized ansatz 




FIG. 7: fcT-factorized results compared to the full CYM cal- 
culation in the dilute limit g^/i^ — 0.2 GeV, m — 0.1 GeV. 



for the gluon multiplicity 

dA^ ~°2 AT 2 



ttR\ N/ - 1 



dyd2kT (27r) 



d^qr 



2Nc 7rg2k^ J (27r)2 
(kT - <lT)'^^fB (kr - Qt 



q|A/;(qT) 



(29) 



where (for nucleus A) 



^AcLt) 



N, 



^ fdhr e*'"^"'^ (Tr[/t(xT +rT)f/(xT) 

(30) 
and similarly for A/"^ with U replaced by V. The tilde 
refers to the adjoint representation, since the adjoint 
correlator is what appears in the A-side of the pA fc^- 
factorization formula. We emphasize that there is no 
derivation of Eq. ([29]) in the "AA"-case, it is an ansatz 
that a) is symmetric in the two nuclei and b) reduces to 
the correct known result in both the "pA" and the "Ap" 
limits, when one of the sources is taken to be weak. 

The formula ([^^ lends itself to simplifications in differ- 
ent limiting cases. In particular, in the large momentum 
limit Ik-rl ^ Qs the integral is dominated by the regions 
q-r ~ and qr ^ k-p. In the first one of these regions 
one can approximate k^ — q^ ~ k^ in the second factor 
and pull it outside the integral; in the second region the 
same can be done to the other factor. This leaves the 
result 



dA^ 



7:R\ N^^ - 1 



dyd2kr (27r) 



2Nr 



T^g 



N^ {\iT)xG^ {x, k|) + A/; (kT)xG^ {x, k^,) (31) 



that is proportional to the dipole cross section or Wilson 
line correlation function M{\<iT) in one nucleus and the 
integrated gluon distribution 



xG{x, krp) 



IkTl ^: 



qr 



(27r) 



q|A/'(qT) 



(32) 
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kp-fact 

kp-fact w/ CO. 1 
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FIG. 8: fcT-factorized results compared to the full CYM calcu- 
lation for the "pA" case g^^.^ = 1.25 GeV, g'^ fXp = 0.008 GeV. 
The dot-dashed curve ("c.o. 1" i.e. "cutoff 1") is the case 
where the momentum from the proton is restricted to be less 
than that of the produced gluon (i.e. nucleus A in Eq. (|29p is 
taken as the proton and there is a cutoff 6 (| kr | — | Qt | ) ) . In the 
dotted curve ("c.o. 2") the cutoff is applied to the momen- 
tum from the nucleus (i.e. B in Eq. H29p is the proton and the 
cutoff is still ^(Ikrl — Iqrl)). The full numerical result shows 
a significant oscillation which is typical for a calculation at a 



finite ending time g fi^r = 60. 



Ny = 20 and m ■ 



0. 



in the other one. At kx S> Qs we have parametrically 
Af{\i.T) ^ Ql/^T (note that A/'(xt) is dimensionless) and 
xG ~ Ql 

Two main variants of the fc^-factorized formula 
Eq. ([29)) have also been used in the literature. 



• The qr-integration is often cut off (see e.g. [28 
^^) with ^dk^l — IqtI); this is asymmetric in the 
two nuclei but does have the advantage of making 
the spectrum IR- finite, which it is otherwise not. 
This cutoff is relatively easy to implement in our 
numerical evaluation, and we shall discuss its effect 
below. 

• Instead of q^A/'(qT) (which we know from the pA- 
case to be the object appearing on the A-side) one 
sometimes replaces the q^ by a l/r|, inside the 
Fourier-transform Eq. pop . This gives an uninte- 
grated gluon distribution <f{q_T) (the "WW" distri- 
bution, see [3,[3,[l3,|3j] ) that is related to the num- 
ber of gluons as defined in LC quantization and be- 
haves like In Iqrl for small momenta. This is closer 
to the KLN ansatz [SJ, [33] '/'(qr) ^ est. /as. Due 
to the logarithmic divergence of the unintegrated 
gluon distribution at small qr this approximation 
is more difficult to treat directly in our numerical 
setup and we do not study it further here. 

Figure [7| shows the comparison of Eq. ([^^[1 with or 
without the cutoff 6'(|kT| — Iqrl) to the full numerical 



kp-fact 
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FIG. 9: ^T-factorized results compared to the full CYM calcu- 
lation in the saturated strong field regime, with g^fi — 2 GeV, 
m = 0.1 GeV Ny — 20 (same configurations as in Fig. [3]). 



result in the dilute regime where one expects agreement. 
One sees that indeed, as can be shown analytically, in 
the dilute case the /cT-factorized expression (without the 
cutoff) is accurate, but the cutoff changes this already 
at quite high momenta. Figure [8| shows the comparison 
between Eq. ([^^ with and without the cutoff and the full 
result in the asymmetric "pA" case. Again we confirm 
the analytical calculation showing that the /cy-factorized 
approximation is good also in the dense-dilute case. Note 
the dependence on whether the momentum that is cut off 
is that of the proton or the nucleus. Figure [9| shows the 
same comparison between Eq. ([29[) and the full CYM re- 
sult for the "AA" case. We see that while the cutoff does 
make the spectrum IR finite, it does so at the expense of 
deviating from the full result already at high momenta, 
kx ^ 2(7^/x. Unlike in the previous dilute cases, also the 
result without the cutoff deviates from the full result for 
kx ^ ff^M- Comparing Figs. [9| and [3| one immediately sees 
that the approximation using the Coulomb gauge fields 
at T = is much more accurate than the fcT-factorized 
one. 



Note that our result on the form of the spectrum in 
/cT-factorization does not invalidate computations where 
/cT-factorization has been used to study the dependence 
of the integrated multiplicity on energy, rapidity, cen- 
trality etc. The integrated multiplicity will still be, by 
dimensional reasons, proportional to Q^. Thus the de- 
termining aspect of these phenomenological applications 
is the dependence of Qs on impact parameter and x, not 
the precise shape of the initial gluon spectrum, which will 
be modified later in the plasma phase. 



VI. CONCLUSION AND PERSPECTIVES 

In conclusion, we have studied numerically the spec- 
trum of gluons in the Glasma fields in the initial stages 
of a heavy ion collision. We have compared the results 
obtained by approximations where one includes the non- 
linear interactions of the gluon fields only in the initial 
condition (at t = 0) to the full numerical GYM calcula- 
tion. The separation between initial and final state ef- 
fects is not a gauge invariant one, and thus our discussion 
naturally involves finding a gauge that minimizes the fi- 
nal state effects. We find that in the Fock-Schwinger -|- 
transverse Coulomb gauge the effect of final state rescat- 
terings on the spectrum is surprisingly small except at 
very small momenta. It would be interesting to see how 
the corrections to the linearized approximation converge 
towards the full CYM result. 

We have also compared our results to those obtained 
assuming /cT-factorization. While in the "pp" and "pA" 
cases the results are the same, as is well known, in the 
fully nonlinear "AA" case /cr-factorization gives a poorer 
description of the gluon spectrum for kx ^ Qs than the 
linear approximation used in this paper, with a marked 
sensitivity to the infrared cutoff. 



CYM 

VU gauge 1 
VU gauge 2 




FIG. 10: Different discretization methods of the VU-gauge, 
see text for the explanation of the labels, g^fj, = 2 GeV, 
Ny = 1 m = 0. Note that the spectrum is multiplied by fcy 
to show the effects at large fcr. 



and fin) can be extracted from Eq. (pS)) . 
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Appendix A: Proton-Nucleus collisions and dilute 
limit 



It was shown in Rcf. [12] that Eq. (P^ evaluated with 
the n, — FC/-gauge fields gives the known fcr-factorized 
formula for the gluon spectrum in the "pA" case. It was 
shown in Ref. [22l that a Coulomb gauge calculation of 
the gluon spectrum in pA collisions gives the same result. 
Let us here briefly show how this comes about evaluating 
our Eq. ((23)) in Coulomb gauge. Assuming nucleus B to 
be a proton one can expand the gauge field to first order 
in p^ , and we get 



A^r^O^ - 


= f^(o)[^l,^L(i)]^!o) (Al) 


-A(i)^=0+ = 


. f^(i)A^ n^^ + n^o)A\ nl^ + n^o^Ai^,^ n}^) 




-^^(i)^'^!o) - ^g^io)d%y 



where fi/Q) = U , 



A' 



(1) 



V'd+'^^' 



(A2) 



(A3) 



^(1) = wf- 



d 



U'^{UA^^,,^U^). (A4) 



We obtain after some algebra 



4'' 



-4(1) 



=0+ 



r=0+ 



-(5^C/)A;(i)C/t-C/4(i)(a'C/t),(A5) 



d'd^ 



^) (uA^.i.^U^) . (A6) 



Plugging Eq. (jASP into (P^ leads to the well known 
/cT-factorization formula for gluon production in proton- 
nucleus collisions. 



Appendix B: Hamiltonian variables used in the 
numerical calculation 

In the numerical calculations it is customary to work in 
a Hamiltonian formalism with the gauge potentials and 
electric fields 



E' 



tAj 



A^ = -r^A" = x+A- - x-A^ 



En 



— A^ri 



(Bl) 
(B2) 

(B3) 



The transverse gauge potential is represented in terms of 
the link matrix 

U^^^-igaA.^ (B4) 

In Refs.[ij, [lal the notation Ay^ = 4> was used, while the 
longitudinal electric field was denoted E^ = tt in [la| and 
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FIG. 11: GYM gluon spectrum at different times. Ny 
20, m = 0.1 GeV (same configurations as in Fig. [S]). 



E^ = p in [13] • In terms of these Hamiltonian variables 
the initial condition Eq. ([5]) for the longitudinal field is 



i-r)\T=0+ 







i?''|.=o+ = -tg[Al,^,Ay 



(B5) 
(B6) 



In terms of the Hamiltonian variables the expression 
for the multiplicity with linearized final state evolution, 
Eq. ([231), reads 



dN 



1 



1 



dyd^kT (27r)2 7rk| 



|kTX^(kT)| 



T = 0+ 



(B7) 



The relation between the notations in Ref. [16| and 
Rcf. [12] is, with ^ on the left and [l2| on the right 
of the equal signs, C/(i) — V^ , C/(2) — U\ Al^. = ^^b 
and Al^-. = —A^^ (the sign is compensated by the op- 
posite sign in the covariant derivative 9^ + igA'f. = 
9^ — igAf^ ^). Between these two references the Wilson 
line is exchanged with its Hermitian conjugate (including 
changing the direction of the path ordering in t he p ath 
ordered exponential); in the conventions of Ref. [ifl| the 
pure gauge field is 



A 



(i)--^^«^'^(i) 



(B8) 



Appendix C: Discretization of Eq. (|23p in the VU 
gauge 

The formulas we need to discretize are Eq. ([M)) : 

A' = V{UA^U^ -A^)V\ (CI) 

and Eq. 1^5^ 



AW [(d'U) A^U^ + UA'^ (d'U^)] V\ (C2) 
which can also be written as 

-4" = eVab (d'Ubc) A'f. 



(C3) 



A straightforward way to discretize Eq. (jCll) is to de- 
fine 



A\^t) 



+ 



V'(xT+JT)(c/(XT+JT)^i(xT)C/t(XT+JT) 
-A^^(xT))f/t(xT+JT) 
F(xT)(t/(xT)A^,(xT)t/t(xT) - A^^(xt))t/^(xt) , 



(C4) 



where Jt is a vector of length a in the j-direction. Here 
the lattice gauge field is really the antihermitian part of 
the link matrix: 

A,,(xt) - -A,\^t) - -^ (UsA^t) - C/I.(xt)) , 

(C5) 
where the link matrices representing the transverse pure 
gauge fields of the individual nuclei are 



Ufi^TJ = C/t(xT)t/(xT + iT) 
Ufi^T) = F^(xT)F(xT + iT). 



(C6) 



In terms of this we then write down the discretized 
derivative needed in Eq. (P5|) as 

e'^d.A' (xt) = e'^ [A' (xT + It) - A^ (xt)] ■ (C7) 

The longitudinal field can is then discretized in the 
same spirit as the transverse one. The two versions 
Eqs. (jC2p and (jC3p lead to different-looking discretiza- 
tions 



A^ = -V{^T + h) 



L/(xt + iT)AB,;(xT)J/^(xT + ir) - U{xt)A^^{xt)U''{xt) 



vH 



■K.T + It) 



+ 2^(^^) 



L/(xT + iT)A^^{yLT)U\-KT + it) - U{yiT)A^^{-KT)U\yiT) 



Ft(xT) (C8) 
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and 



A^ = -V{^t + it) 



[/(xT + iT)A^,(xT)C/UxT) + C/(xT)A^,(xT)[/^(xT + iT)-2C/(xT)A^^(xT)C/t(xT) 



V\yiT + 'iT) 



2^(xt) 



U{^T + iT)A^,(xT)C/t(xT) + ;7(xT)A^,(xT)C/t(xT + It) - 2C/(xt) A^ ,(xT)t/t (xt) 

I 



F^(xt). (C9) 



Of these two we prefer the first one Eq. (|C8p because of 
its relative simpUcity. The combination Eqs. (jC4l) and 
(IC8l) is labeled as "VCZ-gauge 2" in Fig. [TOl 

An alternative discretization method is suggested by 
the observation that the version of the light cone gauge 
proposed in [i2| is equivalent to taking the gauge fields 
in the Fock-Schwinger gauge and performing a gauge ro- 
tation with the product of the Wilson lines of the two 
nuclei. The lattice version of initial conditions for the 
transverse fields is obtained [13, [l^ by solving the link 
matrix Ui{x.T) from the equation 



Tr 



ta ({Ut + C/f ) (l + Uf) - h.c. )] = 0. (CIO) 



Here the pure gauge link matrices U^ ' (see Eq. (|C5p ) 
correspond to the pure gauge fields of the two nuclei sep- 
arately. For SU(2) this equation can be solved in closed 
form, but in the case of SU(3) that we are interested in 
here it is solved numerically by an iterative procedure. 
This link matrix Ui is then used in the initial condition 
for the longitudinal field 



^''(^^) = I^E [(c/.(xT)-l)(^^(xT)-f/^(xT) 



+ {uI{^t-it)-1) (;7,^(xT-iT)-t/f (xt-It))- h.c. 



(Cll) 



The fields obtained from Eqs. (JCIOP and (|C11I) arc then 
gauge transformed with either Q, — UV or ft = VU and 
taking the antihermitian part of the link matrix as the 
gauge field in Eq. ^. F[/-gauge 1 in Fig. \M We 
emphasize that the small difference between these two 
methods only appears at momenta of the order of the 
lattice UV cutoff. We observe that the UV behavior of 
the latter method is closer to the Coulomb gauge curve, 
and it is the one we use in the rest of this paper. 



Appendix D: Multiplicity in the CYM calculation 

In the CYM computations one has to choose some fi- 
nite time T at which to perform the Fourier decompo- 
sition of the fields for computing the spectrum. In the 
boost invariant calculations the interactions of the fields 
become weaker with time, so the result does not depend 
very strongly on t for t 3> l/Qs- The remaining residual 




FIG. 12: CYM gluon spectrum at g^fJ.T — 12 using different 
methods. The solid line labeled "Sum" is the one defined by 
Eq. (jD3|) that we use, unless otherwise stated, in this paper. 
The dashed line ("AE") is the formula (|D2|) that cancels the 
unknown dispersion relation and the dotted one ("EE") the 
one of Eq. ([HB- Ny = 20, m = 0.1 GeV (the field configura- 
tions are the same as in Fig. [3]) . 



dependence is demonstrated in Fig. [TTl which shows the 
resulting gluon spectrum for g^/xr = 4.2, 7.8 and 12. 

Let us then turn to the question of defining the mul- 
tiplicity corresponding to a classical field configuration 
that does not yet evolve completely linearly, as is the 
case in practical numerical computations. The straight- 
forward method used in [16,] is to start from the electric 
field part of the Hamiltonian in Coulomb gauge. The 
electric fields represent approximately half of the total 
energy of the system (which is easily verified in the full 
numerical computation) and because the electric field 
part of the Hamiltonian is quadratic in the canonical mo- 
menta one can obtain a decomposition of the total energy 
of the system into transverse momentum modes. If we 
now assume that each of these modes has a free disper- 
sion relation uj{c[t) = IqtI (or the corresponding lattice 
equivalent in numerical computations), we get our first 
definition of the multiplicity 



dA^ 



1 



1 



dyd'^kr (27r)2 [krl 



+ T7ra(kT)7ra(-kT) 



(Dl) 



The dispersion relation of the interacting theory is. 



12 



however, not free. The imphcations of this observation 
for our situation were observed aheady in the first CYM 
determination of the Glasma multiphcity [ij| . By look- 
ing at correlators of the fields A^ and the corresponding 
electrical fields separately one can numerically determine 



this dispersion relation. It was observed in [l4| that this 
dispersion relation exhibits a mass gap that decreases 
with time as m^ '-^ g^^/r. By taking products of the 
correlators one can then construct a formula for the mul- 
tiplicity where the dispersion relation cancels out: 



J 



dN 



dy d^ky 



f27r 



^ '-ElikT)Eli-kT) 



TTraikT)7Ta{-kT)^TA^ai^T)Al{-kT) + -M^t)M'^t)- (D2) 

I 



The third option, and the one we will use unless oth- 
erwise stated, is based on the formal derivation from a 
reduction formula [36]. Again this is similar in spirit to 



our first definition and amounts to assuming a free disper- 
sion relation uj{\iT) and taking the sum of the canonical 
fields and momenta as 



J 



dN 



dyd^kr 2(27r)2 [w(kT) [t 



--El{kT)Eli~kT)+T7ra(kT)7Tai-kT) 



-aj(kT) 



TAl(kT)A^(- 



-kT) + -0a(kT) 

r 



^ai-'kr) 



(D3) 



When done properly the reduction formula also gives 
an additional (small in practice) contribution that is 
antisymmetric in k^ ^^ — ky. This contribution van- 
ishes when the single inclusive multiplicity is averaged 
over configurations and already for a single configuration 



when the spectrum is averaged over the azimuthal angle 
of k-T- It is, however, irnportant for multigluon correla- 
tions in the Glasma [13, I37h39| . The agreement between 
the three methods is illustrated in Fig. [V2\ 



[1] E. lancu and R. Venugopalan in Quark gluon plasma 
(R. Hwa and X. N. Wang, eds.). World Scientific, 2003. 
[arXiv : hep-ph/0303204 

[2] H. Weigert, Prog. Part. Nucl. Phys. 55 (2005) 461 
_arXi v:hep-ph/0501087 . 

[3] L. D. McLerran and R. Venugopalan, 



[s; 



\Phys. Rev. D49 (1994) 2233 arXi v : hep-ph/9309289] . 
[4] L. D. McLerran and R. Venugopalan, 

\Phys. Rev. D4 9 (1994) 3 352 arXiv:he p-ph/9311205"| . 
[5] L. D. McLerran and R. Venugopalan, 

Thys. Rev. 050^(1994) 2225 arXiv : hep-ph/9402335"] . 
[6] T. Lappi and L. McLerran, 



Nucl. Phys . A772 (2006) 200][arXiv:hep-ph/0602189] . 
[7] D. Kharzeev, Y. V. Kovchegov and K. Tuchin, 

\Phys. Rev 1 )68 (2003) 094013 

arXiv : hep-ph/0307037 . 

J. P. Blaizot, F. Gelis and R. Venugopalan, 

Nucl. Phys. AT43 (2004) T3| [arXiv : hep-ph/0402256 ' . 
[9] F. Gelis and Y . Mehtar- Tani, 

\Phys. Rev. D73 (2006) 0340191 



[10] 

[11] 



arXiv : hep-ph/0512 079] . 

Gelis, A. M. Stasto and R. Venugopalan, 



Eur. Phys. J. C48 (2006) 489 arXiv : hep-ph/0605087 
Y. V. Kovchegov, Nucl. Phys. A692 (2001) 557, 



[12] 

[13] 
[14] 

[15] 

[16] 
[17] 

[18] 
[19] 
[20] 
[21] 



"arXiv : hep-ph/0 0il252| . 
J. -P. B laizot and Y. Mehtar-Tani, 
[ Nucl. Phys. A818 (2009) 97 
arXiv : 0806 . 1422 [hep-ph] . 
A. Krasnitz and R. Venugopalan, 

\Nucl. P hys. B557 (1999) 237 arXiv : hep-ph/9809433"] 
A. Krasnitz and R. Venugopalan, 
\Phys. Rev. Lett. 86 (2001) 1717] 
"arXiv : hep-ph/0007 108 . 
A. Krasnitz, Y. Nara and R. Venugopalan, 



Phys. Rev. Lett. 87 (2001) 192302 
arXiv:hep-ph/0108092 . '^~^~ 



T. Lappi, Phys. Rev. C67 (2003) 054903 

larXiv : hep-ph/0 303076 . 

A. Krasnitz, Y. Nara and R. Venugopalan, 

'Nucl. Phys. A727 (2003) 427" arXiv : hep-ph/0305112] 



T. Lappi, Eur. Phys. J. C5 5 (2008) 285 ; 
arXiv: 071 1.3039 [hep-ph2[^ 
T. Lappi, S. Srednyak and R. Venugopalan, 
JHEP 01X2010) 066| [arXiv : 0911 . 2068 [hep ^^ph]] . 
A. Kovner, L. D. McLerran and H. Weigert, 



Phys. R ev. D52 (1995) 3809 arXiv : hep-ph/9 505320] . 

A. Kovner, L. D. McLerran and H. Weigert, 

Phys. Rev. D52 (1995) 6231, .arXiv:hep-ph/9502289] . 



13 



[22] 
[23] 

[24] 

[25] 
[26] 
[27] 

[28] 

[29] 
[30] 

[31] 



A. Dumitru and L. D. McLerran, 

,Nud. Phys. A700 (2002) 492 arX iv : hep-ph/0105268] . 

P. Romatschke and R. Venugopalan, 

\Phys. Rev. Lett. 96 (200 6) 062302J 

TarXiv:hep-ph/0510121 . 

P. Romatschke and R. Venugopalan, 

Phys. Rev. D74 (2006) 045011 

arXiv:hep-ph/0605045 . 
A. Krasnitz, Y. Nara and R. Venugopalan, 
Phys. Lett. B554 (2003) 2T| [arXi v:hep-ph/ 0204361] . 
A. Krasnitz, Y. Nara and R. Venugopalan, 
Nucl. Phys. A717 (2003) 268 arXiv : hep-ph/0209269| . 
T. Lapp! and R. Venugopalan, 
Phys. Rev. C74 (2006) 05490 51 
' larXiv : nucl-th/060902 1 1 . 

N. Armesto, C. A. Salgado and U. A. Wiedemann, 
\Phys. Rev. Lett. 94 (2005) 022002 

arXiv : hep-ph/04070 18 . 

T. Hirano and Y. Nara, Nucl. Phys. A743 (2004")~305 

arXiv:nucl-th/0404039l. 

H.-J. Drescher, A. Dumitru, A. Hayashigaki and 
Y. Nara, \Phys. Rev. 074^200 6) 044905 
larXiv : nucl-th/0605012l . 
H. J. Drescher and Y. Nara, 



larXiv : nucl-th/ 061 1017 

[32] J. L. AlhaceieTlPh ys. Rev. Lett. 99 (2007) 262301 



[arXiv: 0707 .2545 [hep-ph]] ~ 

[33] J. L. Albacete, N. Armesto, A. Kovner, C. A. Salgado 

and U. A. Wiedemann, 

Thys. Rev. Lett. 92 (2004) 082001J 

arXiv : hep-ph/0307179 . 
[34] D. Kharzeev and M. Nardi, 



\Phys. Lett. B507 (2001) 121 | |arXiv:nuc l-th/001 2025 



Phys. Lett. B523 (2001) 79| 



[35] D. Kharzeev and E. Levin 

'arXiv : nucl-th/0108006 ' 
[36] F. Gelis, T. Lappi and R. Venugopalan, 

Unt. J. Mod. Phys. E16 (2007) 2595J 

[arXiv: 070 8. 0047 [hep-ph] . 
[37] F. Gelis, T. Lappi and R. Venugopalan, 

\Phys. R:ev. D78 (2008) 054020, 

[arXiv: 0807 .13 06 [hep-ph]]7"^ 
[38] F. Gelis, T. Lappi and R. Venugopalan, 

Phys. Rev. D79 (2008) 094017 

■larXiv: 0810. 4829 [hep-ph] . 
[39] K. Dusling, F. Gelis, T. Lappi and R. Venugopalan, 

td. Phys. A836 (2 010) 159| 
Kiv: 091 1.2720 [hep-ph]] . 



\Phys. Rev. C75 (2007) 034905 



